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We present a general definition of damage spreading in a 
pair of models. Using this general framework, one can define 
damage spreading in an objective manner, that does not de- 
pend on the particular dynamic procedure that is being used. 
The formalism is applied to the Domany-Kinzel cellular au- 
tomaton in one dimension; the active phase of this model is 
shown to consist of three sub-phases, characterized by differ- 
ent damage-spreading properties. 
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I. INTRODUCTION 

The concept of damage spreading was introduced in the 
context of biologically motivated dynamical systems by 
Stuart Kauffman |l[ . The question posed is whether the 
phase-space trajectories of two slightly different copies of 
a dynamic system, subjected to the same thermal noise, 
will stay close (or even merge) at long times or, alter- 
natively, will they diverge? Damage spreading made its 
first appearance in the physics literature in the mid eight- 
ies [gr§|, an d attracted considerable interest and atten- 
tion. The main reason behind this initial enthusiasm was 
the hope that damage may spread (indicating chaotic be- 
havior) in some regions of a system's parameter space and 
disappear or heal elsewhere. This possibility intrigued re- 
searchers, since if indeed realized, it would have indicated 
the existence of different dynamic phases in various com- 
plex systems (such as spin-glasses) The initial en- 
thusiasm concerning damage spreading has abated dur- 
ing subsequent years; the main reason being an appar- 
ent lack of an objective, observer-independent measure of 
whether damage does or does not spread in a given sys- 
tem. Even for relatively simple models, such as the two 
dimensional ferromagnetic Ising model, different results 
were obtained when heat bath or Metropolis dynamics 
were used . Both these dynamic procedures are phe- 
nomenological (since they satisfy detailed balance, they 
can be used to generate equilibrium ensembles) and the 
two are equally legitimate to mimic the temporal evo- 
lution of a system in contact with a thermal reservoir. 
If spreading or healing of damage were to indicate some 
intrinsic property of the system, one would not expect 
the result to depend on the details of exactly which phc- 
nomenological procedure was used to generate its dynam- 
ics. 

The purpose of this communication is to pose the 
"right" question; i.e. one which has a well defined ob- 



jective answer. The essence of the argument is to con- 
sider the entire family of dynamic procedures that are 
consistent with the physically dictated constraints of the 
problem. For any particular system one of three possibil- 
ities may hold: 

1. Damage is spreading for every member of the family 
of dynamic procedures 

2. Damage heals for every member of this family 

3. Damage spreads for a subset of the possible dy- 
namic procedures, and heals for the complementing 
subset. 

Hence the only question regarding damage spreading that 
has an unambiguous, observer-independent answer is: to 
which of these three classes a particular system belongs? 

To demonstrate the general concept introduced here 
we studied the simplest dynamic model in which dam- 
age spreading has been observed, the one-dimensional 
Domany-Kinzel (DK) cellular automaton JtJ, for which 
we found the phase diagram presented in Fig. [j]. 
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FIG. 1. Phase diagram of the Domany-Kinzel automaton. 
The active, percolating phase consists of three sub-phases; 
each is numbered according to the damage spreading class to 
which it belongs (see text). 

The DK automaton is a two-parameter model whose 
temporal evolution contains, as special cases, the bond 
and site directed percolation problems. The main point 
made by DK was universality: namely, that the entire 
family of observed transitions of the one-dimensional cel- 
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lular automaton is in the directed percolation universal- 
ity class. DK identified two phases; a "dry" or "frozen" 
phase, in which all initial conditions evolve to the absorb- 
ing state, and an "active" or percolating phase. Some 
years later Martins et al || discovered that in a certain 
region of the active phase damage spreads, and it heals 
elsewhere. More detailed investigations, using simula- 
tions H - [fl2f as well as analytic (mean field) approxi- 
mations fll3|] - jlj] confirmed the existence of this 
"chaotic phase" . Its boundary, however, was shown fllf l 
to depend on the manner in which the dynamic proce- 
dure of the underlying DK model is carried out, while 
the evolution of a single replica is completely insensitive 
to the dynamic procedure. This prompted Grassberger 
|l2f to observe that "it is misleading to speak of different 
phases in the DK automaton... instead these are different 
phases for very specific algorithms for simulating pairs of 
such automata" . This observation is the precise analog of 
the problematic nature of viewing DS as a manifestation 
of a dynamic transition in spin models, where, as men- 
tioned above, it was well known that different dynamics 
that yield identical equilibrium properties can give rise to 
different results for damage spreading. Thus, again, DS 
becomes a "subjective" concept, which is devoid of well 
defined meaning for the DK model, whose phases should 
be determined by the properties of a single evolving sys- 
tem. 

The main purpose of this paper is to point out that if 
one defines the most general family of dynamic rules that 
are consistent with the physics of the problem being stud- 
ied (Sec II), DS has an objective, observer-independent 
meaning. Past work on DS in the DK model is reviewed 
in Sec III and in Sec IV the existence of the three well 
defined distinct phases described in the Introduction is 
established for the DK model by numerical simulations 
and analytical arguments. 

We also tested and confirmed a recent conjecture of 
Grassberger, to the effect that the damage spreading 
transition is in the directed percolation universality class 
]l2fl . Analytical support for this conjecture came so far 
from approximate mean-field arguments fl4|| and an exact 
statement first made by Kohring and Schreckenberg [jl3| , 
who noted that on the p2 = line the dynamics of dam- 
age spreading in the DK automaton is precisely identical 
to the evolution of the DK automaton itself, and hence 
on this line DS is trivially in the DP universality class. 
This being a rather special line, it is of interest to try 
to establish such precise mapping of DS to DP also else- 
where in the pi — P2 plane. In Sec IV A we present such 
an extension. 



II. DAMAGE SPREADING - GENERAL 
FORMALISM 



A. Rules for Legitimate Damage Spreading 
Procedures 

We turn now to present our arguments for the possi- 
bility of defining an observer-independent measurement 
of damage spreading. By this we do not mean that DS is 
reflected in the dynamic behavior of a single system, so 
that Grassberger's observation still holds; DS is a prop- 
erty of a pair of automata.^ It is possible, however, to 
address the lack of objectivity implicit in one's freedom 
to choose the precise algorithm that is used for the evo- 
lution of the pair of replicas. If every observer can pick 
his favorite dynamic rule, get results (on DS) that de- 
pend on the rule used, while no measurement done on 
an evolving single system can differentiate between the 
rules - indeed it appears contradictory to claim that DS 
reflects "phases" of the model that is being investigated. 
Nevertheless such phases can be defined in a precise way. 

To overcome this apparent paradox we formulate quite 
general and physically motivated restrictions on the pos- 
sible dynamic rules that one can use for studying DS. By 
"physical" we mean that the restrictions are dictated by 
the dynamics of the single evolving system. The restric- 
tions are as follows: 

1. The dynamic rules for the evolution of the pair 
of replicas are such that the evolution of a single 
replica is according to its "natural" dynamics. 

2. The transition probability matrix for a site i for 
the pair of replicas can depend only on those sites 
that affect the evolution of site i under the dynamic 
rules of a single system. 

3. The rules that govern evolution for the pair do not 
break any of the symmetries of the single-replica 
dynamics. 

The first restriction simply means that the fact that we 
are watching two systems evolving in parallel should not 
affect the behavior of any one of them. The second con- 
straint means that if the evolution of site i is affected, say, 
only by the states of its nearest neighbors, the relative 
states taken on site i by the two replicas should not feel 
longer range interactions. For example, if site i and all 
its neighbors j are in the same state in the two replicas, 



In this sense DS, as defined for stochastic dynamics, differs 
from dynamics in deterministic nonlinear systems. In that 
case one can find signatures of chaotic behavior in following 
the phase-space trajectory of a single evolving system. Di- 
vergence of two initially neighboring trajectories (indicating 
the existence of a positive Lyapunov exponent) is a compu- 
tationally feasible tool to ascertain the chaotic nature of a 
single system's trajectory, but is not essential to consider two 
replicas in order to define chaos. 



2 



we do not expect damage to be generated at i by a dam- 
aged site which is far away (i.e. not one of the neighbors 
of i). The third rule implies, for example, that if there is 
a left-right symmetry in the evolution of a single system 
the same must hold for the pair of replicas. 

Clearly, the subjectivity in defining the damage spread- 
ing procedure that was described above has now been 
shifted to this point - to selecting the restrictions that 
define which DS procedure is "legitimate" . We do believe 
that there is much less arbitrariness, however, in this kind 
of subjectivity than what was done before - choosing, at 
random, one out of a continuum of physically equivalent 
procedures. 

B. Transfer matrix formalism 

We now introduce a general formalism for damage 
spreading problems in one-dimensional models with par- 
allel updates^. Consider a one-dimensional spin chain 
with L sites and certain boundary conditions to be spec- 
ified. At each site j — 1, ...,L a local spin variable 
(jj = 1, . . . , s is placed which can be in s different states. 
We denote by P t (ai,.. ., ox) the probability to find the 
system at time t in the state {cti, . . . , ox} which is posi- 
tive and normalized by 

E P t (a 1 ,...,a L ) = 1. (1) 

cti,...,o\l 

The time evolution of a system with parallel updates is 
discrete and can be described by a transfer matrix T , 
whose element T^' - '^ is the probability of making 
a transition from state {a%, . . . , ox} at time t to state 
{n, . . . ,Tl} at time t + 1. The corresponding discrete 
time master equation has the form 

p t+1 (n,...,r L ) = E T l\';::^ L Pt^i,---^L) 

a u ...,a L 

(2) 

The conservation of probability (|l|) implies that 

E T ^'.v:,^ = i (3) 

n,...,Tt 

for all configurations {<7i, . . . , ox}, which means that the 
elements in each column of the transfer matrix add up to 
one. 

In what follows we assume that each site is updated 
simultaneously 2 and that an update at some site j de- 
pends locally on n sites of the previous configuration, 



2 This is suitable for various cellular automata. For systems 
that satisfy detailed balance one must not update interact- 
ing spins simultaneously. Generalization of our formalism to 
random sequential update is straightforward. 



say {aj- m , . . . , <Tj +m } where m = ^yl- This means that 
the transfer matrix factorizes: 

L 

T Tl,...,T L _ TT rpTj / 4 \ 
t7l,...,<Tl J_J_ M tJj- m ,...,<7j + m Y*l 

3=1 

Each factor Tjj_ mj ... !(Tj . +m conserves probability sepa- 
rately: 

s 

E T «i- m ,~,°i +m ( 5 ) 

We now consider two replicas S and S' of the same sys- 
tem. Denote by {<7i, . . . , e>x} and {ti, . . . , 7x} the states 
of S and by {a[, . . . , a' L } and {t[, . . . , t' l } those of S' . In 
order to define simultaneous temporal evolution of these 
systems, one has to generalize the transfer matrix: 

rpTl,...,TZ, ^ IJ\T1,...,TL M ,T 1 ,...,T L 

cri,...,ax o"i,...,o\£, ; tx\ ,...,(tC 

rpTj _^ rp T J ; T j 

CTj_ m ,...,(Tj +m o-j_ m ,...,(Tj + m ; a' j m ,...,a' j+m 

The transfer matrix for the total system (S, S') is re- 
stricted by the requirement that each replica should 
evolve as before, i.e. integrating out the degrees of free- 
dom of one of the replicas results in the previous transfer 
matrix of the other replica: 

s 

VT Tj:Tj ' , , =Tl , (6) 

/ j Cn'_ 7Tl ,...,<Tn'_i_ Tr i : (7 ,a. i a a . , v / 

T ]= \ 
S 

V T Tj ; Tj ' , = T T i (7) 

t'. = 1 

3 

These restrictions already imply probability conservation 
for the total system: 

s 

y T Ty ' T ' ] , , =i (8) 

A * J III I , J^IIL , j_ m l 5 j + m 

Tt ; t' — 1 

In order to study damage spreading, a further restriction 
is imposed: Once both replicas reach the same state (no 
damage), their temporal evolution is identical: 

rp T i • T 3 _ rp Tj tt 3 ; r' (q\ 

J-a-j— m ,...,a-j +m ;a-j- m ,...,aj +m — ± aj- m ,...,a j + m u \ v ) 

We define damage as the Hamming distance between the 
states {(7i, . . . , ox} and {crj , . . . , a' L } at any given time, 
that is the fraction of sites for which Oj ^ o'y 

A = zE A . = zEa-W (10) 

3=1 3=1 

The restrictions 
@; (0); an( i @ impose dependences among the s 2n+2 
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matrix elements T 3 ' 3 , , . The num- 

crj_ m ,...,crj +m ; a } _ m ,—,cr j+m 

ber of independent degrees of freedom can be counted 
as follows. First notice that because of eq. (||) only 
2s — 1 of the 2s equations in (Q) and (]?]) are independent. 
Thus, even though for any given initial configuration 
{a-j-m, . . . , (Tj +m ; Oj_ m , . . . , cr' j+rn } the number of differ- 
ent final configurations is s 2 , the number of independent 
matrix elements is s 2 — (2s — 1) = (s — l) 2 . On the other 
hand eq. M) implies that for s n of s 2n possible initial con- 



figurations {<jj 



5 0~ j+m ■ 



j+m 



} the 



trix elements are already defined. Furthermore the whole 
system is symmetric under exchange of the replicas which 
gives another factor |. Thus the total number of inde- 
pendent degrees of freedom of the transfer matrix is 



1 



„2n 



a") (a- l) 2 



(11) 



C. Algorithmic implementation 

In a numerical simulation the temporal evolution de- 
scribed by the transfer matrix T can be realized as fol- 
lows. At each time step all sites are updated indepen- 
dently, i.e. <jj is replaced by a new value tj according 
to probabilistic rules which depend only on the previ- 
ous configuration <Tj- m , . . . , <7j+ m . Since in this paper 
we discuss only two-state models, let us from now on re- 
strict our attention to the case s = 2, with <7j = 0, 1 (the 
generalization to s > 2 is straightforward). 

We introduce a stochastic binary 
variable r a . 



, , that denotes the value assigned to 
site j in one update, given the state of its neighborhood 



' j-\-m ■ 



(12) 



The numbers r a . mi ... ta -. +m are generated probabilisti- 
cally in some procedure, such that updates at different 
sites or different times are uncorrelated. Furthermore, 
their expectation value, averaged over many realizations 
of random numbers, is given by the corresponding matrix 
element of the transfer matrix: 



(13) 



Usually this procedure is implemented by generating a 
random number z from a uniform distribution in the in- 
terval < z < 1, and comparing it with the transition 

r — 1 

probability T CT J _ mi ... i(T - +m , to get the assigned value of 



i if • 

otherwise 



(14) 



This prescription is not unique, however, and there are 
many other possibilities to select J-<T,-_ m ,...,<T,- +m in a man- 
ner that satisfies eq. (|l3|). This constraint specifies 



the temporal evolution of a single system. The exis- 
tence of correlations between any two random variables 
r a , and ?v a i will play no role in the 

evolution of a single system, since only one of the two 
will be used for any given update. The situation is dif- 
ferent, however, in damage spreading problems in which 
the two different replicas of the system are evolving si- 
multaneously. The updated state of the pair of replicas 
is governed by the joint transition probabilities 



Prob(r 3 



,)=T T 



,(7 j + m ; CT. 



(15) 



Therefore correlations between the two random variables 
r a , r a i will influence the temporal evolution of the com- 
bined system and, therewith, the properties of dam- 
age spreading. Such correlations are contained in the 
two-point correlation functions of the random variables, 
which are given by the elements of the combined transfer 
matrix for both replicas: 



T 



r 3 -=l;r<=l 



(16) 



In this formalism the restriction (||)-(0), stating that each 
replica separately evolves in the same way as the original 
system, is satisfied automatically. The same applies to 
the second restriction (^|) which ensures that in case of 
'no damage' both systems evolve in parallel. 

Obviously, the number of degrees of freedom specifying 
damage spreading is just the number of two-point corre- 
lation functions. In the case of two-state models there 
are i2™(2™ — 1) such correlations, which agrees with the 
number of degrees of freedom counted in eq. ( pi] ) . 

Three-point correlations do not affect the evolution 
of a pair of replicas, because in each update only two 
random variables, r a ._ CT ., and r ff < CT / , are 
used. However, they would start to play a role in dam- 
age spreading problems with three replicas. Generally, 
fc-point correlations will be felt in systems consisting of 
at least k replicas. 



III. DAMAGE SPREADING IN THE DK MODEL: 
A BRIEF REVIEW 

In this section we review briefly past work on dam- 
age spreading in the DK automaton. We emphasize the 
manner in which DS was calculated by various authors, 
and the manner in which different ways of defining DS 
are embedded in the general formal framework of Sec II. 

The DK automaton is defined as follows: a binary vari- 
able (Ji(t) = 0, 1 characterizes the state of site i at (dis- 
crete) time t. a — 1 means that the site is wet or active, 
whereas a = means that it is dry. The automaton 
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evolves by a parallel update rule, which can be stated, 
using the notation of Sec II, as follows: 

T Ti = T Ti (17) 

!0 if (Tj_i = a i+ i = 
Pi if <n-i <r»i (18) 
p 2 if <7j_l = CT i+1 = 1 

That is, the state of site i at time t + 1 depends only on 
the states of its two neighbors at time t; only wet sites 
can give rise to a wet site, with probabilities p\ if one 
neighbor was wet and P2 if both were wet. 

Using the notations introduced in Sec II the transi- 
tion probabilities in the DK model are defined by the 
one-point expectation values of three stochastic binary 
variables 

(nil) = (no) =pi , (m)=P2- (19) 

This model has a dry phase and a wet phase, separated 
by a transition line which has been determined with high 
accuracy by various numerical methods. In spite of its 
simplicity, the model has not been solved exactly, except 
for the special line pi = 1 , [pi , |l7j . At all points on 
the phase boundary, except the special line, the transition 
to the active or wet phase is characterized by directed 
percolation (DP) exponents. 

Damage spreading properties between two replicas are 
controlled by correlations between the random variables: 

(roi m) = (r w r n ) = & 

(rot no) = (20) 
(mi r w ru) = 7 

According to the arguments discussed in Sec II, only one- 
and two-point functions enter the transfer matrix, which 
means that 7 is an irrelevant parameter in the present 
problem. 

Martins et al || were the first to address the issue of 
damage spreading in the DK model. Two nearly identical 
initial configurations were allowed to evolve on two repli- 
cas, using the same random numbers for both (the precise 
meaning of this statement will be explained below). They 
discovered that the active phase contains in fact two re- 
gions; one in which damage spreads and its complement, 
where it does not. The boundary between these regions 
was subsequently determined with increasing accuracy by 
Zebende and Penna ||, by Martins et al p0|] , Rieger et 
al and Grassberger Independently, mean- field 

type approximations of varying complexity were also used 
to study the DS problem pl|Jl^ -[l5[| . An interesting ob- 
servation, first made by Kohring and Schreckenberg [[l3[ , 
was to the effect that the position of the "phase bound- 
ary" between the DS and non-spreading regimes depends 
on the manner in which damage is generated. The orig- 
inal scheme of Martins et al used a single uniformly dis- 
tributed random number < z < 1 for the two replicas: 
using the definitions of Sec II this means that the choice 



''01 = r 10 = 8(p 1 - z) , rn = 8{p 2 - z) (21) 

was made, which can be expressed as 
(3 = (romo)= Pi 

a=<n B iru)=Min(pi,p 3 ) (22) 

The dynamical process is generated by setting 

o-j(t + 1) = n = r CTs _ 1(Ti+1 
a[{t + l) = T l l =r <i<+i 

For reasons that will be evident shortly, we call this dy- 
namic process, that uses a single random number, maxi- 
mally correlated. 

Kohring and Schreckenberg recognized the fact that 
one could, in principle, use two different random numbers 
to determine and r 4 ', if at least one of the two neighbor 
sites was damaged at time t. In fact they studied DS 
using two different random numbers zoi and zn, their 
DS procedure has two fully correlated binary variables 
(roi and r\o) and two uncorrelated ones (roi and rn): 

roi = no = 6{pi - zoi) , m = 0(p2-zu) (23) 

the correlations being 

(3 = (roi no)= Pi 

a = (r i ni)= (r i)(m) = Pi P2 (24) 

The dynamics generated by using on the first replica Ti = 
ro- i _ 1 a i+1 and t[ = r a > on the second gave rise to a 

shift of the original "phase boundary" (as obtained with a 
single random number, eq. ((2^))- As discussed in Sec II, 
the evolution of a single replica is completely insensitive 
to whether one or two random variables are used in the 
dynamic procedure, which prompted Grassberger p^ ] to 
make his observation quoted in the Introduction. 

Finally we note that Grassberger has formulated re- 
cently |12j a conjecture, which is a natural extension 
of previous statements H,E9|,0 regarding universality 
of directed percolation transitions for models with non- 
symmetric absorbing states p0|| . According to this con- 
jecture damage-spreading transitions should be in the 
universality class of directed percolation fl^| , provided 
some general conditions are satisfied. The DK model is 
a natural candidate to test this conjecture because of its 
simplicity, ease to simulate and our precise knowledge of 
the existence of a DS transition and its location. Grass- 
berger presented numerical evidence for his conjecture, 
which we also confirmed and extended. We also show 
below that in a region of the p\ , P2 plane one can map 
DS exactly to the DK model and hence onto DP. This re- 
sult is an extension of a statement first made by Kohring 
and Schreckenberg ]13] ], who noted that such a mapping 
holds on the P2 = line. 
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IV. TRUE PHASES IN THE DK MODEL 

As discussed in Sec III, the most general dynamic rule 
that can be defined for two replicas of the DK automaton, 
in accordance with these constraints, has two degrees of 
freedom or parameters, a and (3. As it turns out (see 
Appendix |A|), the possible values that a and /3 can take 
are restricted by requiring that all transition rates in the 
transfer matrix have to be positive. For any value of 
Pi,p 2j the range of allowed values of the parameters a 
and P is given by 

max(0,pi +p 2 - 1) < a < mm(p 1 ,p 2 ) 
max(0, 2a - p 2 , 2p\ — 1, 2pi — 1 — 25 < P < Pi 

(25) 

There are three important special cases, namely those of 

• maximal correlations: a = min(pl,p2), (3 = p\ 

• no correlations: a = Pip 2l P = P\ 

• minimal correlations: a =pi + p 2 — 1, 

In the case of minimal correlations, the values listed 
above hold only in the region 2pi + p 2 > 2 (see Ap- 
pendix |a|) . 

A. Exact results 

We turn now to show that for p 2 /2 < p\ < \~p 2 /2 the 
damage spreading process can be mapped exactly onto a 
directed percolation process. Kohring and Schreckenberg 
p3| have shown that such a mapping holds on the line 
p 2 = 0. Clearly, their choice of parameters (^) is a par- 
ticular case of our damage spreading procedure, which 
is the most general one that satisfies rules 1-3 listed 
above. Therefore we find a wider (two-dimensional) re- 
gion in the p\ , p 2 plane in which such a mapping is pos- 
sible. To see this, let Aj = 1 — 5 CTii<T < be the damage at 
site i. By Po(Aj = 1 | <7j_i<jj_|_i; cr-_ 1 cr- +1 ) we denote the 
probability to generate a damaged site for a given initial 
configuration in a particular update. These probabilities 
are listed in Table I, in which we introduced for brevity 
the notation 

X = 2pi - 2/3 Y = pi + P2 - 2a (26) 

In general, the probability for generating damage on site 
i depends on the previous states of both replicas, i.e. on 
((Tj_i(Tj_|_i; Oi_iOi + 1); knowledge of Aj_i and Aj+i does 
not suffice to determine A, at the next time step. Thus 
damage spreading itself cannot be seen as an independent 
process. We may, however, pose the following question: 
under which conditions will damage spread as if it were 
generated by an independent process? That is, when do 
we have 



P D (Ai\o- i -io- i+1 ;a' i _ 1 a' i+1 ) = Pd(A,- | A,-_i, A, + i) 

(27) 

In order to satisfy this condition, any two entries in 
Table I, that correspond to the same initial damage 
{Aj_iAi + i}, should be equal. For example all four initial 
configurations 

= ({11; 00}, {10; 01}, {01; 10}, {00; 11}) 

have the same initial damage {Aj_i,Aj+i} = {1,1}- In 
order to satisfy eq. (p7|) , the four entries (p 2 ,X,X,p 2 ) 
must have the same value, i.e. we must have p 2 = X. A 
similar consideration leads to the condition Y = pi ; that 
is, we must have 

Pd(1|00) = 

Pd(1|01) = Pd(1|10)=,Pi=Y (28) 
P D (1\U) = f>2 = X 

Note that these are precisely the update rules of the DK 
process. Using the definitions (^6|), we see that the cor- 
relations must satisfy 

a= y> / 3 = Pi-y- ( 29 ) 

Since the correlation parameters are restricted by 
eq. ( p5| ) , the allowed range for p\ and p 2 in which these 
conditions can hold is a triangle in the phase diagram: 

P2/2 < pi < l- P2 /2 (30) 

To summarize: we have proved that within this trian- 
gle we can find correlations a, (3 such that the damage 
spreading process follows the dynamical rules of a single 
DK automaton. 
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TABLE I. Probabilities P D (A.i = l\a i -i(Ti + i;a' i _ x a' i+x ) 
for the generation of damage in the DK model. X and Y are 
defined in eq. (26) 
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Say we have a line in the (pi,P2) plane that lies within 
this region. For every point (Pi,P2) on this line we can 
find a, (3 values for which DS evolves precisely like a DK 
automaton with parameters Since part of the 

transition line of the DK model (from dry to wet phase) 
lies in the triangle (|30|), on any trajectory that crosses 
this part of the phase boundary we will observe a damage 
spreading transition precisely at the DP transition and 
with DP exponents (provided we chose <S, (3 according to 
eq. ((2^).) In particular, this holds for the line P2 = 0, 
as discovered in [Ol ; note that for P2 = their choice of 
correlations, eq. fl24l ) precisely satisfy eq. (p9|) . 

B. Results from comparing probabilities 

Other useful results can be obtained by comparing 
probability tables of different pairs of automata. The ba- 
sic idea is that by increasing (decreasing) all probabilities 
in Table I, damage spreading will be more (less) likely. 
More precisely, if a pair of DK automata described by 
parameters p*,p*,, a*, ,9* exhibits damage spreading, we 
expect that any other pair of automata with parameters 
Pi,P2,a,f3 satisfying 

Pi < Pi 

P 2 < P2 

Pi+P*2~ 2a* < pi + p 2 - 2a (31) 
2p\ - 2(3* < 2pi - 2(3 

exhibits damage spreading as well. Vice versa, if damage 
heals in a pair of automata described by pl,p*,, a*, (3*, 
then for any other pair with pi, P2, a, (3 obeying 

Pi > Pi 

P2 > P2 

pl+pl-2a* >pi+p 2 -2a (32) 
2p\ - 2(3* > 2p x - 2(3 

we expect damage to heal. 

Although these statements are very plausible, we were 
not able to prove them rigorously. However, we per- 
formed various numerical tests which turned out to be 
consistent with the inequalities stated above. 

Because of these inequalities the boundaries between 
the three regions in the phase diagram correspond to ex- 
tremal correlations a and (3. For example, if at a point 
{P11P2) damage spreads in a model with maximal corre- 
lations a max = mm(pi , P2) and (3 max = p lt then eq. (|3l| ) 
implies that damage spreads also for every a and (3 in 
the allowed range (p5|). This, however, means that the 
point (pi,j?2) belongs to region 1 in the phase diagram. 
Therefore the phase boundary of region 1 coincides with 
the DS transition line for maximal correlations. Similarly 
one can use eq. ([}2|) to show that the phase boundary 
between regions 2 and 3 coincides with the DS transition 



line for minimal correlations. It turns out (see Fig. g) 
that this line lies entirely in the region 2p\+p2 > 2 so that 
minimal correlations are well defined (see Appendix [^) . 

Alternatively one can compare the probabilities for 
generating damage in a pair of DK automata to the prob- 
abilities of generating a wet site in a single DK automa- 
ton. To this end one simply has to use the same in- 
equalities setting a* = p*./2 and (3* = p\ - j>*./2. For 
example, if p\ and p*. represent a point in the wet phase 
of the DK phase diagram, then for all pairs of automata 
parametrized by pi,P2, ct, (3 and satisfying 

Pi < Pi 

P2 < P2 

P*i<Pi+P2-2a (33) 
p*2 < 2pi - 2/3 

damage will spread. On the other hand, if p\ and p* 2 
belong to the dry phase of the DK model then in all 
pairs of automata with 

Pi > Pi 

P*2 > P2 

p$>Pi+p 2 -2a (34) 
p*2 > 2 V i - 2(3 

damage does not spread. 

As an illustration of eq. ( p3| ) consider the point Mi in 
Fig. |: Setting p\ = p\ w 0.809 and p* 2 = we obtain 
the conditions 

Pi < Pi , Pi < Pi +P2 - 2a , < 2pi - 2f3 

(35) 

Using the bounds (p5| ) we find from these inequalities 
that in the triangle pi ~ P2 > p\ damage spreads with 
certainty for any a in the allowed range. In Fig. I this 
triangle is indicated as a shaded region. 

C. Terminal points of the phase boundaries 

We turn now to derive, using the arguments introduced 
above, a few exact results concerning the phase bound- 
aries for minimal and maximal correlations. As explained 
above, these boundaries are the transition lines between 
the damage-spreading phases 1,2 and 3 described in the 
Introduction and shown on |l|. In order to make our ar- 
guments easier to follow, we present in Fig. || all the lines 
and special points that are mentioned. 

Let us consider first the case of maximal correlations. 
That is, for every point (^1,^2) in the phase diagram 
we assign a — min(pi,p2) and (3 = p\ and look for 
the boundary B max between the region in which dam- 
age spreads and the one in which it doesn't. Denote the 
wet-to-dry transition line of the DK model by B we t (see 
Fig. |). 
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We now prove that B max and B we t can intersect only 
at the point Mi = (pi,0), where B wet intersects the pi 
axis. On this axis maximal correlations correspond to 
the choice a = P2 = and /3 = p\ , which also satisfy the 
conditions (29), i.e. the damage spreading process can 
be mapped onto a single DK model. Therefore we know 
that a DS transition will occur precisely at p\ = pi, so 
that at all points (pi > p1,P2 — 0) we must have DS. 
On the other hand, as we will now show, for maximal 
correlations there cannot be DS on any point on B we t, to 
see this, note that (for the region of interest, pi > 1/2) 
maximal correlations imply the inequalities 



the segment of B we t that lies above the intersection point 
is in the region where pi +p 2 /2 > 1 and this immediately 
leads (for minimal correlations) to the inequalities 



Pi + P2 - 2a = \pi -p 2 \ <pi 
2 Pl - 2/3 = < p 2 



(36) 



P 2 




FIG. 2. Schematic phase diagram, displaying the special 
lines and points that are referred to in the text. 

For a point {pi,p 2 ) on B we ± these inequalities, when 
used together with eq. (|34|) , imply that damage does not 
spread. Furthermore, since ([36|) are (for p2 > 0) strict 
inequalities, B wet lies inside the no-spread region. The 
point Mi is on the boundary between this region and a 
region (containing the p\ > pi axis) in which damage 
does spread; Therefore Mi must lie on B max . 

We turn now to the phase boundary B m in for minimal 
correlations, and show that it terminates at the point M2 , 
where B wet intersects the line C, given by p\ +P2/2 = 1. 
First notice that on C minimal correlations correspond 
to a = 1 — pi and {3 — 2p\ — 1. Therefore (see eq. (|29|)) 
the DS process is equivalent, on C, to a DK model, so 
that as we move on C, keeping minimal correlations, a 
DS transition occurs at the point M% = (p*,P2)i where C 
intersects B we t- As before, we show next that at all points 
on B we t with P2 > p\ damage doesn't spread; hence B m i n 
must pass through M 2 . To prove the last claim note that 



Pi + P2 - 2a = 2 - pi - p 2 < pi 
2 Pl -20 = 2- 2 Pl < p 2 



(37) 



According to eqs. (p4|), this implies that on the DK tran- 
sition line in the region p\ + P2/2 > 1 damage does not 
spread. Therefore M2 is the terminal point of B m in- 

Having located the endpoints M\ and Mb, we now turn 
to the opposite end of the lines Bmax and B m i n . Note that 
for p\ = 1 the bounds ( p5| ) collapse to a = P2,j3 = 1, 
i.e. maximal and minimal correlations are identical and 
hence B max and B m i n meet at some point M3 on the 
Pi — 1 line. The three special "multicritical" points dis- 
cussed above determine the topology of the phase dia- 
gram for DS. In order to obtain high precision quantita- 
tive information about the location of the transition lines 
we performed numerical studies of damage spreading in 
the DK model. 



D. Numerical results 

In order to obtain accurate numerical estimates for the 
critical parameters of models with absorbing states one 
usually has to let the system evolve for extremely long 
times Q. Grassberger overcame the difficulty posed by 
long transients and obtained good statistics by simulating 
n replicas of the same system in parallel, using simple 
bit manipulations on computers with unsigned words of 
length n [ fL2| . Using this multi-spin encoding method he 
measured the decay in damage on a one- million site chain, 
allowing it to evolve for hundreds of thousands of time 
steps. Because of the improved statistics he was able to 
determine the critical exponents for damage spreading at 
a particular transition point with high accuracy. 

Another method to determine the critical point effi- 
ciently is the so-called gradient method which was intro- 
duced by Zebende and Penna || . In this method a gradi- 
ent in pi and P2 is arranged along the chain. The values 
of the parameters at the two end-points of the chain are 
chosen to be in different phases, i.e. on different sides 
of the transition point. This allows the critical point to 
be determined by measuring the average location of the 
boundary of the active (damaged) cluster. 

In the present work we used a combination of multi- 
spin encoding and the gradient method. In combining 
these methods, a number of problems emerged which we 
solved as follows: 

1. In order to measure the damage spreading transi- 
tion point, one has to find the first position (ap- 
proaching from the non-spreading phase) where 
damage occurs. Simulating n = 64 lattices in par- 
allel, this has to be done for each of the 64-63/2 = 
2016 pairs of replicas. To do this one has to set 
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up a 64 x 64 table in order to keep track of dam- 
aged pairs. Moreover, one has to scan the words 
bit by bit which makes it impossible to use parallel 
bit manipulations. The large amount of CPU time 
needed for this process usually kills the advantage 
one gains from the multi-spin encoding. In order to 
solve this problem, we used a simplified search al- 
gorithm which is based on fast bit operations. The 
price we pay is that only w 75% of all possible 
pairs are taken into account^. We proved that the 
error of this method docs not bias the measurement 
of the transition point. 

2. Zebende and Penna started each run with a sin- 
gle damaged seed located somewhere on the chain. 
It is not clear whether the choice of the location 
influences the results. In order to circumvent this 
problem, we used initial conditions with randomly 
distributed damage all over the chain. 

3. The gradient method is a finite-size simulation and 
therefore boundary conditions may play an im- 
portant role. In the work of Zebende and Penna 
the boundary conditions can be understood as dry 
walls and it is not clear to what extent they affect 
the measurements. In order to minimize this ef- 
fect, we created, on a chain of 2N sites, a gradient 
with reflection symmetry (pi(i) — p\{2N — i) and 
Pi(i) = Pi{2N — i)) and measured the boundary 
of the damaged cluster on both sides. We expect 
finite-size effects to be less important for these pe- 
riodic boundary conditions. 

The phase diagram of the DK automaton, obtained us- 
ing the multiple lattice gradient method, is presented in 
Fig. [I]. First, we verified numerically the prediction that 
larger correlations correspond to smaller damage and vice 
versa. This was done by scanning the (a, (3) space for 
various points in the (pi,P2) plane. Next we determined 
the DS transition lines for minimal and maximal corre- 
lations. Typical gradient values of 1.2 • 1CP 5 were used 
for lattice sizes L — 8192 and upwards. A transient pe- 
riod of at least 2L was followed by an averaging period 
of L time steps. For (pi,P2) near the the transition lines 
longer transient times were used. The terminal points 
of the phase boundaries were determined with high ac- 
curacy. Using a chain with L = 16384 sites, gradients 
down to 1.22 • 10" 6 and transients of 231072, we mea- 
sured the following critical values at these special points: 
pi = 0.8087(5) (on the p 2 = line); p% = 0.3130(5) (on 



the pi = line). The new triple point was located at 
p\ = 0.744(10), p% = 0.526(10). This was done without 
using our analytic result that identified this point as the 
crossing of C with B we t; the value of p\+P2/2 = 1.007(15) 
agrees with the predicted value (of 1) for points on C. 

Measuring the density of damage along the gradient of 
the chain, we could estimate the density exponent (3. At 
the terminal points we found (3 — 0.302(30) for P2 = 
and (3 = 0.296(30) on the p\ = 1 line. We also measured 
the exponent at a point {p\,p2) = (0.85,0.35) which lies 
inside phase 3. This was done by crossing the DS phase 
boundary while varying the correlations a and (3, yielding 
the value (3 — 0.279(10). All results are in fair agreement 
with the expected density exponent of directed percola- 
tion (3 = 0.277(1) ||, ||. 

V. SUMMARY 

We have rules that a most general damage spreading 
procedure should satisfy. These rules are most natu- 
ral: they ensure that the evolution of a single replica 
is not affected by the fact that two replicas are evolving 
simultaneously; that the range of damage spreading does 
not exceed the range of interactions in the original sin- 
gle model and that the two evolving replicas respect the 
symmetries of the model. These rules can be cast in a for- 
mal setting, using transfer matrices. Using these formal 
definitions we were able to parametrize the most gen- 
eral damage spreading procedure for any given model in 
terms of correlation coefficients between various stochas- 
tic binary variables. Thus we are considering all possible 
damage spreading procedures and identify different dam- 
age spreading phases in terms of the manner in which 
this complete set of procedures behaves. Three possible 
phases can occur; one in which damage spreads for all 
allowed procedures, one in which it does not spread for 
any procedure and the third, in which for some proce- 
dures damage spreads while for others it does not. 

These ideas were implemented for the Domany-Kinzel 
automaton, for which the three phases were identified, 
using a combination of numerical and analytic methods. 
We have shown that in an extended region of the model's 
parameter space damage spreading can be mapped onto 
the evolution of the DK automaton itself. This observa- 
tion supports Grassberger's recent conjecture to the ef- 
fect that damage spreading is in the directed percolation 
universality class. This was also confirmed by numerical 
tests (performed in regions where the above mentioned 
mapping does not hold). 



3 In this approximation, a given replica is declared to be dam- 
aged at site i if the majority of the other replicas is in a dif- 
ferent state. Instead of keeping track of damaged pairs of 
replicas, we recognize only single replicas from where damage 
originated. This amounts in dropping statistically 25% of all 
possible pairs. 



APPENDIX A: GENERATION OF CORRELATED 
RANDOM VARIABLES 

In this Appendix we explain in detail how correlated 
random variables roi, r\Q and rn, that govern the evolu- 







tion of the DK-model, can be generated. We also prove 
the allowed ranges for a and (3, given in eq. (^5j) . Finally, 
we explain the manner in which minimal correlations are 
given by the expression presented in Sec. V. 

Since in each update 7*01, rio, and ru can be either 
zero or one, there are eight possible combinations. By 
^r m ,r uu rn wc denote the (positive) probability to gener- 
ate the combination {roi, rio, ru}- These probabilities 
are normalized by: 

T000 + 7T001 + 7T010 + 7T011 + 7T100 + VT101 + 7Tno + 7T1H = 1 

(Al) 

Once the probabilities T rol . ri0!ril are known, the num- 
bers roi, riQ, and ru can be generated by taking one 
uniformly distributed random number < z < 1 and se- 
lecting one of the eight possible outcomes according to 
the probabilities T roiiri0jril . The three-point correlation 
functions can be represented in terms of the 7r's: 

Tin = (roiriorn) 

ttiio = (roirw (1 - ru)) 

Appropriate combinations of these probabilities are re- 
lated to two- and one-point functions and therewith to 
the entries of the transfer matrix; for example, 



ttiio + tiu = (roirio) = Tq{ 
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(A2) 



Collecting all identities of this type, we obtain seven 
equations: 

TTlOO + 7T101 + THO + 7Tin = P\ 
ToiO + T011 + ""110 + Till = Pi 
T001 + Toil + T101 + Tin = P2 

Ton + Tin = a (A3) 
Tioi + Tm = a 

TllO + Tm = P 
Till = 7 

Together with the normalization ( |Al[ ) these equations 
determine all probabilities T roliri0jril : 

tooo = 1 - 1 + P + 25 - 2pi - p 2 

tooi = P2 - 2a + 7 
toio = tioo =_pi-a-/3 + 7 (A4) 
Ton = Tioi = a — 7 

TllO = P~ 7 

Till = 7 

Since all n have to be positive, we obtain six inequalities: 

max(0, 2a - p 2 ) < 7 < a (A5) 
max(7, 7 - 1 - 2a + 2p\ + p 2 ) < P < 7 + Pi ~ a (A6) 

For a given choice of the parameters p\ , p 2 these inequal- 
ities imply restrictions on the correlation parameters a, 



and 7. The allowed range of these parameters can be 
derived as follows. First let us consider the restrictions 
on a. Eq. ( |A5| ) implies that < a < p 2 whereas eq. 
( |A6| ) leads to the condition pi + p 2 — 1 < a < pi . Both 
of them can be combined by requiring 



max(0,pi + P2 — 1) < a < min(pi,p2) 



(A7) 



For a given a in this i nterv a l th e maximal ranges of (3 
and 7 are given in eqs. ( A5)-( |A6| ). However, since we do 
not explicitly use the three-point correlation parameter 
7, we are only interested in the maximal range of j3. This 
range can be o bta ined by inserting the extremal values 
for 7 into eq. (A6), that is 7 = max(0, 2a — pi) on the 
l.h.s. and 7 = aon the r.h.s. Thus for a given a in the 
range (A7) the corresponding maximal range of (3 is: 



max(0, 2a - p 2 , 2p x - 1, 2pi +p 2 - 1 - 2a) < (3 < pi 

(A8) 



In other words, if a and (3 satisfy eqs. (|A7|)-(A8), we are 
able to find some 7 such that all probabilities T roliri0jril 
are positive. 



APPENDIX B: MAXIMAL AND MINIMAL 
CORRELATIONS 

In order to determine the phase boundaries B max and 
fimin, one has to identify situations where the correla- 
tions are extremal. For given parameters pi,p 2 maximal 
correlations (minimal damage) simply corr espon d to tak- 
ing the upper bounds of the intervals (A7) and (A8): 



a = min(pi,p 2 ) 



P = Pi- 



(Bl) 



In case of minimal correlations (maximal damage) the 
situation is more complicated since th e co rrelation pa- 
rameter a appears on the l.h.s. of eq. (A8) with a neg- 
ative sign. Therefore by increasing a, the minimal value 
of (3 may decrease which makes it impossible to predict 
where damage is maximal. However, we can show that 
in the triangular region above line C in Fig. |^, where 
2pi +p 2 — 2 > 0, this problem does not arise. In this re- 
gion eq. (A7) implies that a > p\ +p 2 — 1 and there with 
2a > 2pi + 2p 2 — 2 > p 2 . Hence according to eq. ( [A5| ) 
the minimal value that 7 can have is 2a — p 2 . Inserting 
this value into eq. (A6) we get the inequality 



max(2ci — p 2 , 2p\ — 1) < [3 < pi 



(B2) 



which replaces eq. (A8) in the specified triangle. In this 
inequality a occurs with a positive sign on the l.h.s. and 
therefore the case of minimal correlations is well defined: 



Pi + Pi - 1 , 



P = 2 Pl - 1 , 



(B3) 
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